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ABSTRACT 

Context. In the presence of strong density stratification, turbulence can lead to a large-scale instability of a horizontal magnetic 
field if its strength is in a suitable range (within a few percent of the turbulent equipartition value). This instability is related to a 
suppression of the turbulent pressure so that the turbulence contribution to the mean magnetic pressure becomes negative. This results 
in the excitation of a negative effective magnetic pressure instability (NEMPI). This instability has so far only been studied for an 
imposed magnetic field. 

Aims. We want to know how NEMPI works when the mean magnetic field is generated self-consistently by an a 2 dynamo, whether 
it is affected by global spherical geometry, and whether it can influence the properties of the dynamo itself. 

Methods. We adopt the mean-field approach which has previously been shown to provide a realistic description of NEMPI in 
direct numerical simulations. We assume axisymmetry and solve the mean-field equations with the PENCIL CODE for an adiabatic 
stratification at a total density contrast in the radial direction of s=s 4 orders of magnitude. 

Results. NEMPI is found to work when the dynamo-generated field is about 4% of the equipartition value, which is achieved through 
strong a quenching. This instability is excited in the top 5% of the outer radius provided the density contrast across this top layer is at 
least 10. NEMPI is found to occur at lower latitudes when the mean magnetic field is stronger. For weaker fields, NEMPI can make 
the dynamo oscillatory with poleward migration. 

Conclusions. NEMPI is a viable mechanism for producing magnetic flux concentrations in a strongly stratified spherical shell in 
which a magnetic field is generated by a strongly quenched a effect dynamo. 
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1. Introduction 

The magnetic field of stars with outer convection zones, includ- 
ing that of the Sun, is believed to be gene rated by differential ro- 



by dy namo action in turbulent convection (IGuerrero & K apvla. 

1DJ. 



Moffa till978tlParkei 



tation a nd cyclonic cony e ction (s ee, e.g 

19791 IZeldovich et ail Il983t [Brandenburg 

& Subramanian 

2005). The latter leads to an a effect, which refers to 



an important new term in the averaged (mean-field) induc- 
tion equation, quantifying the component of the mean elec- 
tromotive force that is aligned with the mean magnetic 



field (see, e.g.. ISteenbeck et all 119661 I krause & RadTeR Il980t 



IBrandenburg et all 120131) . However, what is actually observed 
are sunspots and active regions, and the description of such 
phenomena is not par t of conventional mean-field dynamo the- 
ory (see, e.g.. iPriestl 119821 IStixL 11989): Ib ssendriiverl l2003t 
ICallv et all 120031: IStenflo & Kosovichevl 120121) ! 



Flux tube mode ls 
Spiegel & Weissl [19801: 



1984 



Dikpati & Charbonneau, 



dParkeri 119551 119821 
Spruitl I198U ISchtissler et al 
19991) have been used to explain the 



1994 



formation of active regions and sunspots in an ad hoc manner. 
It is then simply assumed that a sunspot emerges when the 
magnetic field of the dynamo exceeds a certain threshold just 
above the bottom of the convection zon e and for the duration 
of about a month (Chatt eriee et alll2004l) . Such models assume 
the existence of strong magnetic flux tubes at the base of the 
convection zone. They require magnetic fields with a strength o f 
about 10 5 gauss dD' Silva & Choudhurii [1991 lArlt et all 120051) . 
However, such strong magnetic fields are difficult to produce 



Another possible mechanism for producing magnetic 
flux concentrations is the negative effective magnetic pres- 
sure instability (NEMPI), which can occur in the presence 
of strong density stratification, i.e., preferentially near the 
stellar surface, on scales encompassing those of many tur- 
bulent eddies. D i rect n u merical sim u lation s (DNS; see 
IBrandenburg et all 1201 U iKemeletall 1201 2al). mean-field 
simulations (MFS; see Brandenburg, Kleeorin, & Rogachevskii, 
| 2010fc IB randenbu rg et ail 12012b IKemeletall 120 1 2b ; 
Kapyla et al., 2012) , and ea rlier analytic studies ( Kleeorin et al. , 
1989, 1990 . 119961 iKleeorin & Rogachevskiit 11994 : 
iRogachevskii & KleeorinL [2007) now provide conclusive 
evidence about the physical reality of NEMPI. However, there 
remain open questions to be answered before it can be applied 
to detailed models of active regions and sunspots formation. 

In the present paper we take a first step toward combining 
NEMPI, which is well described using mean-field theory, with 
the a effect in mean-field dynamos. To study the dependence 
of NEMPI on the magnetic field strength, we assume that a is 
quenched. This allows us to change the magnetic field strength 
by changing the quenching parameter. We employ spherical ge- 
ometry using axisymmetry, i.e., azimuthal derivatives are ne- 
glected. Furthermore, a is a pseudo-scalar that changes sign 
about the equator, so we assume that a is proportional to cos 8, 
where 9 is colatitude (Roberts! Il972l We arrange the quench- 
ing of a such that the resulting mean magnetic field is of ap- 
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propriate strength to allow NEMPI to work. This means that the 
effective (mean-field) magnetic pressure has locally a negative 
derivative wi t h respe ct to increasing normalized field strength 
dKemel et all l2012b). so the mean toroidal magnetic field must 
be less than about 20% of the equipartition field strength. 

The choice of using spherical geometry is taken because the 
dynamo-generated magnetic field depends critically on the ge- 
ometry. Therefore, to have a more realistic field structure, we felt 
it profitable to carry out our investigations in spherical geome- 
try. Guided by the insights obtained from such studies, it will in 
future be easier to design simpler Cartesian models to address 
specific questions regarding the interaction between NEMPI and 
the dynamo instability. 

In the calculations presented below we use the PENCIL 
which was us ed in DNS of m agneto-hydrodynamics in 
spherical coordinates (iMitra et al.ll2009l) and also in earlier DNS 
and MFS of NEMPI. Unlike most of the earlier calculations, we 
adopt an adiabatic equation of state. This results in a stratifica- 
tion such that the temperature declines approximately linearly 
toward the surface, so the scale height becomes shorter and the 
stratification stronger toward the top layers. This is done to have 
a clear segregation between the dynamo in the bulk and NEMPI 
near the surface, where the stratification is strong enough for 
NEMPI to operate. The gravitational potential corresponds to 
that of a point mass. This is justified because the mass in the 
convection zone is negligible compared to that below. The goal 
of the present work is to produce reference cases in spherical ge- 
ometry to assess the dependence on numerical resolution, and to 
look for new effects that are due to spherical geometry. We begin 
by describing first the basic model. 



2. The model 

The evolution equations for mean velocity U, mean density p, 
and mean vector potential A, are 



d~A 
~dt_ 
_DU 
1 ~Dt 



UxB + aB-ri T J, (1) 
Jx B + V(q p B 2 /2ti )-is T Q-pVH, (2) 
~ (3) 



where D/Dt = d/dt + U ■ V is the advective derivative, 
H = /i+$ is the mean reduced enthalpy with h = c p T being the 
mean enthalpy, T oc p 7 " 1 the mean temperature, 7 = c p /c v is 
the ratio of specific heats at constant pressure and constant den- 
sity, respectively, $ is the gravitational potential, ~p is the mean 
density, rfr = T)t + V ar, d v t — v t + v are the sums of turbulent 
and microphysical values of magnetic diffusivity and kinematic 
viscosities, respectively, 

-Q = p (V 2 C7 + i VV • U + 2SV lap) (4) 

is a term appearing in the viscous force and Sy = h{Uij + 
Ujj) — g&y V • U is the traceless rate of strain tensor of the 
mean flow, J = V x -B/Mo is the mean current density, /zq is 
the vacuum permeability, and the second term, V(g p -B 2 /2^o), 
on the right-hand side of Equation © determines the turbu- 
lent contribution to the mean Lorentz force. Here, q p depends 



on the local field strength (see below). Note that this term en- 
ters with a plus sign, so a positive q p corresponds to a sup- 
pression of the total turbulent pressure. The net effect of the 
mean field leads to an effective mean magnetic pressure p e g — 
(l — q p )B 2 /2/iq, which becomes negative for q p > 1, which can 
indee d be the case for magnetic Reynolds numbers well above 
unity ((Brandenburg et al., 20 1^). 

Following IKemel et alJ d2012ch . the function q p (J3) is ap- 
proximated by: 



1 



+ 070! pi + p 2 



(5) 



where q p o, /3 P , and /3* = (3 p q p0 are constants, (3 = |J5|/_B Gq is 
the modulus of the normalized mean magnetic field, and B cq — 
y/HoP Urms is the equipartition field strength. 

For NEMPI to work, there has to be a mean horizon- 
tal magnetic field whose strength is in a suitable range such 
that the derivative, dp c s/df3 2 , is negative somewhere within 
the computational domain. Unlik e the Cartesian cases investi- 
gated in earlier work ([Brandenburg. Kleeorin. & Rogachevskiil 
l2010tlBrandenburg et al 1 120 121: IKemel et all 1201 2^ 7 where it is 
straightforward to impose a magnetic field in a sphere, it is eas- 
ier to generate a magnetic field by a mean-field dynamo. In that 
case, we have to include a term of the form aB in the expres- 
sion for the mean electromotive force [see the second term on the 
right hand side of Equation ([TJ]. When the mean magnetic field 
is generated by a dynamo, the resulting magnetic field strength 
depends on the nonlinear suppression of the dynamo. We assume 
here a simple quenching function for the a effect, i.e., 



Q?o cos 9 

1 + Oa/3 2 ' 



(6) 



where Q a is a quenching parameter which determines the typical 

— 1/2 

field strength to be of the order of Q a B eq . The value of Q a 
must be chosen large enough so that the nonlinear equilibration 
of the dynamo process results in that dp e ff / d-B is indeed nega- 
tive within the computational domain. In analogy with the f3 p pa- 
rameter in Equation © we can define a parameter f3 a = Q a '', 
which will be quoted occasionally. 

The strength of the dynamo is also determined by the dy- 
namo number, 



C a = a R/rfi 



(7) 



For our geometry with 0.7 < r/R < 1, the critical value of 
C a for the onset of dynamo action is around 18. The excitation 
conditions for dipolar and quadmpolar parities are rather close 
together. This is because the magnetic field is strongest at high 
latitudes, so the hemispheric coupling is weak. In the following 
we restrict ourselves to solutions with dipolar parity. We adopt 
the value C a — 30, so the dynamo is nearly twice supercritical. 

As mentioned before, our gravitational potential $ is that of 
a point mass. We define $ such that it vanishes at a radius r*, 
i.e., 



$(r) = ~GM I - - — 
r r* 



(8) 



http://pencil-code.googlecode.com 



The radial component of the gravitational acceleration is then 
g = —GM/r 2 . We adopt an initially adiabatic stratification with 
c p T = — $(r), so T vanishes at r = r*, whose value has to be 
chosen some distance above r = R. For our reference model 
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Fig. 2. Dependence of B I1US (solid lines) and U rms (dashed 
lines) on time in units of t]t / R 2 for q p o = (black), 5 (blue), 
10 (red), 20 (orange), 40 (yellow), and 100 (upper black line for 
B rms ). The results for U lms depend only slightly on q p o, and this 



1.00 only when the dynamo is saturated. 



Assume that this equation also applies to the current case where 



Fig. 1. Initial stratification of density and inverse scale height 
for r±/R = 1.001 (strongest stratification), 1.01, 1.05, and 1.1. 
The dashed lines mark the position of the reference radius r rc f = 



H p depends on r, and setting k = H p0 
rate is 



XH, 



0.95 R, where p/p « 0.0068 and H p (r) 
The dotted line marks the value of r\ t //3*u r 



p0 



H 



pO 



H p0 by definition. ^ Uims Hp ^ Uims H, 



ms R-po • 



pQ 



the normalized growth 



(12) 



we use r+/R = 1.001. The radius r* determines the density 
contrast. The pressure scale height is given by 



H P (r) 



r(l - r/rv) 



1 



(9) 



In Figure[T]we compare therefore H p o/H p with r) t / f3*u Tms H P Q 
and see that the former exceeds the latter in our reference model 
with r+/R = 1.001. This suggests that NEMPI should be ex- 
cited in the outer layers. 

For the magnetic field, we adopt perfect conductor bound- 
ary conditions, on the inner and outer radii, ro = 0.7 R and R, 
respectively, i.e., 



dA r 



= A s = Aa, = 0, on r = r , R. 



where n = 1/(7—1) = 3/2 is the poly tropic index for an 
adiabatic stratification with 7 = 5/3. The density scale height 

is H p = r(l - r/r*)/n. Table □ gives the density contrast for ° n the P ole and the ec l uator ' we assume 



different values of r*. The initial density profile is given by 

p/p = (-^/nc%) n . (10) 

Radial profiles ofp/po are shown in Figure [T]for r+/R varying 
between 1.1 and 1.001. 

The analytic estimate of the growth rate of NEMPI, A, based 
on an isothermal lay er with H p = H p = const is given by 
dKemel et all l2012b) 



dA r 



As 



dA^ 
89 



0, on (9 = 0° and 90°. 



(13) 



(14) 



H„ 



Vtk 2 



(ID 



Table 1. Dependence of the density contrast on the value of r*. 
Note that R — r = 35 Mm corresponds to r/R = 0.95. 



*/R 



H p (top) /R 



1.100 3.6 x 10 
1.010 4.0 x 10 
1.001 4.0 x 10 



HpQ I R Pmax / Pmin 

0.052 TlTZTni 
0.023 
0.019 



1.4 x 10 1 
2.9 x 10 2 
8.9 x 10 3 



Since our simulations are axisymmetric, the magnetic field is 
conveniently represented via and A&. In particular, contours 
of r sin OA^, give the magnetic field lines of the poloidal mag- 
netic field, Bpoi = V x (A^cf)). The magnetic field strength 
is given either in units of the local equipartition value, B eq (r), 
or in units of the equipartition value B cq o = B cq (r IC [) at the 
reference radius r rc f = 0.95 R. 

In all cases presented in this paper, we adopt a numerical res- 
olution of 256 x 1024 mesh points in the r and 9 directions. This 
is significantly higher than what has been used previously, even 
in mean field calculation s with stratification and h ydrodynamical 
feedback included; see Brandenbu rg et al.1 (Tl992). where a reso- 
lution of just 41x81 meshpoints was used routinely. In principle, 
lower resolutions are possible, but in some cases we found cer- 
tain properties of the solutions to be sensitive to the resolution. 

3. Results 

In our model, the dynamo growth rate is about 170 t]t/R 2 - 
Although both dynamo and NEMPI are linear instabilities, this 
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Fig. 3. Meridional cross-sections of B$/ B cq (color coded) to- 
gether with magnetic field lines of B po \ for different stratifica- 
tion parameters r* and Q a = 10 3 . The dashed lines indicate the 
latitudes 70.3°, 73.4°, 75.6°, and 76.4°. 



is no longer the case in our coupled system, because NEMPI de- 
pends on the magnetic field strength, and only in the nonlinear 
regime of the dynamo does the field reach values large enough 
for NEMPI to overcome turbulent magnetic diffusion. This is 
shown in Figure [2] where we plot the growth of the magnetic 
field and compare with runs with different values of q p o. For 
q p o = 100 we find a growth rate of about 270 t/t/R 2 , which we 
associate with that of NEMPI. 

Let us now discuss the resulting magnetic field structure. We 
begin by discussing the effects of varying the stratification. To 
see the effect of NEMPI more clearly, we consider a somewhat 
optimistic set of parameters describing NEMPI, namely q p o — 



100 and (3 p = 0.05, which yields /3* = 0.5; see Equation (0. 
This is larger than the values 0.23 and 0.33 found from numeri- 
cal simulat ions with and without sm all-scale dynamo action, re- 
spectively (Brandenburg et aill2012l) . The effect of lowering the 
value of q p o can be seen in Figure [2] and will also be discussed 
below. We choose Q a =^1000 for the a quenching parameter 
so that the local value of B^ / B cq near the surface is between 10 
and 20 percent, whi ch is suitable for the excitation of NEMPI 
(Kernel et all [2012b). Meridional cross sections B^/B^q to- 
gether with magnetic field lines of B po \ are shown in Figure [3] 
Note that a magnetic flux concentration develops near the sur- 
face at latitudes between 70° and 76° for weak and strong strat- 
ification, respectively. Structure formation from NEMPI occurs 
in the top 5% by radius, and the flux concentration is most pro- 
nounced when r* < 1.01. 

Next, if we increase the magnetic field strength by making 
Q a smaller, we see that the magnetic flux concentrations move 
toward lower latitudes down to about 49° for Q a = 100; see 
Figure |4] However, while this is potentially interesting for the 
Sun, where sunspots are known to occur preferentially at low 
latitudes, the magnetic flux concentrations become also weaker 
at the same time, making this feature less interesting from an 
astrophysical point of view. For comparison with the parameter 

/3 p = 0.05 in Equation ((5) we note that /3 a = Q a takes the 
values 0.1, 0.07, 0.04, and 0.03 for Q a = 100, 200, 500, and 
1000, respectively. Thus, for these models the quenchings of the 
non-diffusive turbulence effects in the momentum and induction 
equations are similar. 

Also, if we decrease q p o to more realistic values, we expect 
the magnetic flux concentrations to become weaker. This is in- 
deed borne out by the simulations; see Figure[5] where we show 
meridional cross-sections for q p o in the range 40 < q p o < 100 
for Q a = 10 3 . This corresponds to the range 0.32 < /3* < 0.5. 

For magnetic weaker fields, i.e., for larger values of the 
quenching parameter Q a , we find that NEMPI has a modifying 
effect on the dynamo_in thaHt can now become oscillatory. A 
butterfly diagram of B r and B^ is shown in Figure|6] Meridional 
cross-sections of the magnetic field at different times covering 
half a magnetic cycle are shown in Figure Q It turns out that, 
at sufficiently weak magnetic field strengths, NEMPI causes the 
oscillatory solutions with poleward migrating flux belts. The rea- 
son for this is not well understood, but it is reminiscent of the 
poleward migration observed in the presence of weak rotation 
dLosada et al.L 120121) . Had this migration been equatorward, it 
might have been tempting to associate this with the equatorward 
migration of the sunspot belts in the Sun. 

Finally, we discuss the change of kinetic, magnetic, and cur- 
rent helicities due to NEMPI. We do this by using a model 
that is close to our reference model with r+/R = 1.001 and 
Q a = 1000, except that q p o = in the beginning, and then at 
time to we change it to q p o = 100. The two inverse length scales 
based on magnetic and current helicities, 



ku = 



L, A ■ BdV\ LJ-BdV 

Jv r _ and k c =Mo ^ — , (15) 



increase by 25%, while the inverse length scale based on the 
kinetic helicity, 



f v W ■ UdV 
IyU 2 dV ' 



(16) 



drops to very small values after introducing NEMPI, see e.g. 
Figure [8] Here, W = V x U is the mean vorticity. This be- 



4 



S. Jabbari et al.: Surface flux concentrations and spherical o? dynamo 





Fig. 4. Meridional cross-sections for different values of Q a , for 
r* = 1.001. The dashed lines indicate the latitudes 49°, 61.5°, 
75.6°, and 76.4°. 



Fig. 5. Meridional cross-sections for different q V Q parameters in 
the range 40 < g p o < 100 for Q a = 10 3 . The dashed lines 
indicate the latitudes 68°, 72.5°, 75.7°, and 76.3°. 



havior of is surprising, but it seems to be associated with an 
increase in kinetic energy. The reason for the increase of the two 
magnetic length scales, on the other hand, might be understand- 
able as a consequence of increasing gradients associated with the 
resulting flux concentrations. 

4. Conclusions 

The present investigations have shown that NEMPI can occur 
together with the dynamo, i.e., both instabilities can work at the 
same time and can even modify each other. It was already clear 
from earlier papers that NEMPI can only work in a limited range 
of magnetic field strengths. Hence we have adopted a simple a 



quenching prescription to arrange the field strength in the desired 
range. Furthermore, unlike much of the earlier work on NEMPI, 
we have used an adiabatic stratification instead of an isothermal 
one. This implies that the pressure scale hight is no longer con- 
stant and now much shorter in the upper layers than in the bulk 
of the domain. This favors the appearance of NEMPI in the up- 
per layers, because the growth rate is inversely proportional to 
the pressure scale height. 

There are two lines of future extensions of the present model. 
On the one hand, it is important to study in more detail the in- 
terplay between NEMPI and the dynamo instability. This is best 
done in the framework of a local Cartesian model which is more 
easily amenable to analytic treatment. Another important exten- 



5 



S. Jabbari et al.: Surface flux concentrations and spherical o? dynamo 




Fig. 6. Butterfly diagram of B r (upper panel) and (lower 
panel) for Q a = 10 4 , r* = 1.001, u = 11.3 r) t / R 2 . 

sion would be the inclusion of differential rotation. At the level 
of a dynamically self-consistent model, where the flow speed 
is a solution of the momentum equation, diffe rential rotation 
is bes t implemented by including the A effect ( Riid igerl Il980i 
1989). This is a parameterization of the Reynolds stress that is 
in some ways analogous to the parameterization of the electro- 
motive force via the a effect. Mean-field models with both a 
and A effects have been considered before (Brandenbur g et all 
119921) . so the main difference would be the additional parameter- 
ization of magnetic effects in the Reynolds stress giving rise to 
NEMPI. In both cases, our models would be amenable to verifi- 
cation using DNS by driving turbulence through a forcing func- 
tion. In the case of a spherical shell, this can easily be done in 
wedge geometry where the polar regions are excluded. However, 
in that case the mean-fie ld dynamo solution s are oscillatory with 
equatorward migration (Mitr a et aU 1201 Ob . At an earlier phase 
of the present investigations we have studied NEMPI in the cor- 
responding mean-field models and found that NEMPI can revert 
the propagation of the dynamo wave from equatorward to pole- 
ward. However, owing to time dependence, the effects of NEMPI 
are then harder to study, which is why we have refrained from 
studying such models in further detail. 
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